{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import pylab as pl\n",
    "import itertools\n",
    "pl.rcParams['figure.figsize'] = 4,4\n",
    "\n",
    "from gmm_base import *\n",
    "from simulation import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7efcb9113630>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQQAAAD4CAYAAAAKL5jcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2deXxV1bX4vzvzPJGBhCQkYQhDmGdEpVYtzq21vzqgdrRqta3v2VZ9fbbP5/P1te9j+1qt1loVa0UrtdVaioJDHRiUQSCAYRJCSIAQIIQMZNq/P/a54d7LHc7NPffck2R/P594p3XPWh7uWWfttddeW0gp0Wg0GoCYaBug0Wicg3YIGo2mD+0QNBpNH9ohaDSaPrRD0Gg0fcRFS3Fubq4sKyuLlnqNZsiyYcOGo1LKPF+fRc0hlJWVsX79+mipdyz1J9oBKMpKDksmkjbZrT8UnGybUxBC7Pf3WdQcgsY3d734MQAvfmteWDKRtMlu/aHgZNsGAtohOIw7LxhjiYyVeOuzW38oONm2gYCIVqXizJkzpR4yaDT2I4TYIKWc6eszHSE4jNqmNgBKh6WEJRNJm+zWH4iuri7q6uro6OgAoLunF4C4WD2BlpSURHFxMfHx8aa/ox2Cw/j+ss1A4DGwGZlI2mS3/kDU1dWRnp5OWVkZQgj2NJ4CYFReWpQtiy5SSpqamqirq6O8vNz097RDcBh3XTTWEhkr8dZnt/5AdHR09DkDgIKMpChb5AyEEAwbNozGxsaQvqcdgsOYWzHMEhkr8dZnt/5guJwBQFqi/km7cD8vZtEDLYexp/FUX9gbjoyVeOs7tGYpn9bW2qY/FDq6eujo6om2GQMW7RAcxn0vb+W+l7eGLWMlHvo6Wxn++q2sfukXtukPhYMn2jloFCcNRhYuXBjRgj4dXzmMHyyqtETGSjz09XYDcGGJrSaYZniAHMLj/9zD5OJM5o/K7Xtv9Z6jbKlr5tbzR9lhnuPREYLDmDEyhxkjc8KWsRJf+gpiW2zTHwqpiXGk+skjTC7O5I7nN7F6z1FAOYM7nt/E5OLMsHS2trZy2WWXMWXKFKqqqnjxxRd54IEHmDVrFlVVVdxyyy246n0WLlzIXXfdxXnnncf48eP56KOPuPrqqxkzZgw/+tGPANi3bx/jxo3j5ptvZvLkyVxzzTW0tbWdpfeNN95g3rx5TJ8+nS996UucOhX+MFI7BIdRc6iFmkOBLzYzMlbioc/4YbceO2Sb/lAIlEOYPyqXR66fxh3Pb+LhN2q44/lNPHL9NI+IoT+sWLGCoqIiNm/eTHV1NYsWLeKOO+7go48+orq6mvb2dl577bU++YSEBN59911uvfVWrrrqKh599FGqq6t55plnaGpqAqCmpoZbbrmFLVu2kJGRwW9+8xsPnUePHuXBBx9k1apVbNy4kZkzZ/Lwww+H9f8B2iE4jvtfqeb+V6rDlrESX/qOHj5gm/5QCJZDmD8ql8VzSvnVW7tZPKc0bGcAMGnSJFatWsUPf/hD3nvvPTIzM3n77beZM2cOkyZN4q233mLbtm198ldeeWXf9yZOnEhhYSGJiYlUVFRw4IA6ryUlJZxzzjkALF68mPfff99D59q1a9m+fTvnnHMOU6dOZcmSJezf73fNkml0DsFh3HfpeEtkrMRTn4oQiuJbbbXBLIWZgesQVu85ynPravnOBaN5bl0tc0cNC9spjB07lg0bNrB8+XLuvfdeLr74Yh599FHWr19PSUkJP/nJT/oqKQESExMBiImJ6Xvuet3drXI03lOG3q+llFx00UUsXbo0LNu90RGCw5hSksWUkqywZazEl774jibo7bXNBrOkJMSRkuD7PufKGTxy/TT+5eLKvuGDK6fQX+rr60lJSWHx4sXcfffdbNy4EYDc3FxOnTrFsmXLQj5mbW0ta9asAWDp0qUsWLDA4/O5c+fywQcfsHv3bgDa2trYuXNnWP8foCMEx7GtvhmAiUX+E11mZCJmk2sxnOyF9uOQ6qwipfZOdYdN9uEUttQ1e+QMXDmFLXXNYUUJW7du5fvf/z4xMTHEx8fz2GOP8de//pVJkyZRVlbGrFmzQj7m+PHjWbJkCd/61rcYM2YMt912m8fneXl5PPPMM1x33XWcPn0agAcffJCxY8OrItWrHR3Gl3+r7gqB1gmYkYmYTW3H4GdGbfztayHf3uGLNzt27GD8+DM2DIa1DPv27ePyyy+nujr8PJH3+QG92nFAcf8VEyyRsRK/+lobgeg6BG+KguQQNIHRDsFhmBkG2DVU8KnPPaJsDW3hjB34GioMNMrKyiyJDvqDTio6jM0HTrD5wImwZazEr77W8JJxkaCts5s2I4+gCZ2gDkEIUSKEeFsIsUMIsU0I8V0fMkII8SshxG4hxBYhxPTImDv4eWj5Dh5aviNsGSvx1OfsCKGhuYOG5o7gghqfmImvuoF/lVJuFEKkAxuEECullNvdZC4Bxhh/c4DHjEdNiDxwVZUlMlbiV9+pI7baYYYRuttyWAR1CFLKBqDBeN4ihNgBjADcHcJVwLNSTVmsFUJkCSEKje9qQqByeLolMlbioc8jh+C8IUNSfGy0TRjQhJRDEEKUAdOAdV4fjQDca1nrjPe8v3+LEGK9EGJ9qJ1chgob9h9jw/5jYctYiac+Zw8ZWk9303pa5xD6i2mHIIRIA/4MfE9KedL7Yx9fOavAQUr5hJRyppRyZl6ez41jhjw/W1HDz1bUhC1jJX71OdAhHDrZwaGTOofQX0zN0Qgh4lHO4I9Sypd9iNQB7ivki4H68M0bejx09SRLZKzEQ58xZOiNTSTGgQ4hWjmErVu3cvnll3PPPfecVVU4kDAzyyCA3wM7pJT+1le+CtxkzDbMBZp1/qB/jMpLC1plZ0bGSnzpi0krgM5T0Hn2Ov1okhQfG5U8wqRJk3jhhRd49tlnbddtJWYihHOAG4GtQoiPjffuA0oBpJSPA8uBS4HdQBvwVetNHRqs3avWwwdqZGpGJnI2qQjhVHw2adRC21FIKLXFDjOcMvIH0Wi2mp+f77HMeSBiZpbhfXznCNxlJPBtq4wayvxipVqxFmidghmZiNlkDBlqWpKYASqPkOUQh/CPe4g9YNyzrIoShk+CS35qSvSee+7h9OnT7N+/n5EjR1qj32YGfp3nIOPn10yxRMZKfOmrHDUKtq9zXC1CYlx0im9XrFjR10pt27Zt2iForMHM9mh2b6HmqU9FCGm5xepli4NaqV3y06jU4nd0dPCDH/yAV199laeffprq6mouvfTSKFgSPnotg8N4f9dR3t8VuODHjIyV+NK3uz1VPTl12DY7zNDS0UVLR5etOh988EFuuukmysrKmDRpUtQWJlmBjhAcxq/f2gXAgjH+G3aYkYmYTUYO4c2dxxidkuusCAE40qKahaQnmd/gNBxqampYuXIlH3zwAaBmGx566CFbdEcC7RAcxi++PNUSGSvxpe/Ls0ph63DHRQgl2fYOpyorK1m3bp3Ha1cLtYGIdggOo8hEYY0ZGSvx1KcihKyUBEgrcFyEkBClpOJgQZ89h/FOzRHeqQmcuTcjYyUe+owhwyeHWiB9uOMcQjRyCIMJHSE4jMfe2QPAwsr8sGQibdPbNY2Mm1QArUdU9+WY6N1bpJR9bcrtziE4mf70S9UOwWH8+vpplshYiac+9SO7YU4pxHarvR7bmiAtOovVkpKSaGpqYtiwYQghKM2xN4fgVKSUNDU1kZQUWo9J7RAcRn568H9AMzJW4ktfRnICJBql06cORc0hFBcXU1dXh15OfzZJSUkUFxeH9B3tEBzGqu0qa3/hhIKwZCJmkxGGbqtvYeKkMUqg5bAq8Y0C8fHxlJeX9722+9wMNrRDcBi/e28vEPgHbUYm0jb9c+cRJs43Xrc4Z2Gr3edmsKEdgsN4bPEMS2SsxFOfihBuml8OacPVW6ecM9Ng97kZbGiH4DByUhMskbESD33GkCEtMQ7ikyApUw0ZHILd52awoesQHMaK6gZWVAcOwc3IWIkvfVsOqv0eSS90VIRg97kZbOgIwWE8/cE+ABZVFfqXef/ToDJW4mmTihDe23mUyZdhVCs6J0Iwc/40/tEOwWH87mafe3B68Hzv9+mq+n/A/MgbhG+bvnZuhXqSPhz2r7HFDjOYOX8a/2iH4DAyTFTYxZ7YT+yRzTZYo/CwycghJMcbP520AjVkkBJEwMZatmDm/Gn8o3MIDuNvm+v52+bADau7e3toathniz3g26aNtcfVk/Th0NMJ7cdtsycQZs6fxj86QnAYz63dD8AVU4r8ynT1SHpO1Nllkk+bPtjTxHRQSUVQtQgpObbZ5A8z50/jH+0QHMYzX50dVCYpLoaknmO2LSrysMkYMnzrvFHqdaZRGtt8EAomRtyWYJg5fxr/aIfgMJITgncLFkjo7VIdj9MjX5Hny6aEOOO9DONOfPJgxO0wg5nzp/GPziE4jL9squMvmwIPB7p7jWWtNl2EnjYp3R/tV3s1kDYcRIxjHIKZ86fxj3YIDuOFDw/wwocHAsr09PSqJyftSZ75smn1HmPz19g45RRssiUYZs6fxj96yOAwnvvGnKAyCXEx0I1tF6GHTUYO4dufGXPmvcwR0OyMu7KZ86fxj3YIDiM+NnjQJrB3yOBpk9IdF+NWc5BRBIe322JLMMycP41/9NlzGC+tP8BL64MMGfpyCPZECL5sWvPpsTMvMoqVc+pHyy6rMXP+NP7RDsFhLNtQx7INgcPvnl57cwgeNhkX/bq97g6hCLraoOOELfYEwsz50/hHDxkchpkNXBNiBfRg25DBl03fu3DsmReZI9Rj80FIzrbFJn/YtQHuYEVHCAOZk/VRCNN96MswHIJDZho0/Uc7BIex9MNaln5YG1Cmp7eX7thk6DkNbccCylpuk+GA3t/TdEagzyFEvxbBzPnT+Ec7BIfx2pZ6XtsS+E4rJRwRRsdjGy5CXzZt2O+WL0grcExxkpnzp/GPziE4jD9+Y25QmbgYKCoZBZ/WqTC9cLKNNqkI4bufdatDiI1Ti5wcMGQwc/40/tERwkAls0Q9Nkdpis2790FGkWOKkzT9RzsEh/GHNfv4w5p9AWV6paS6OQFiE21xCB42GTmEf9Z4bYySMcIREYKZ86fxj3YIDmPVjiOs2hFkI1cp2dvUrpYen4h8As2XTZtdTVZdZBarCCHKxUmmzp/GLzqH4DCWfC34ev4YAVdOGQEHS21xCJ42qQv+O58d6ymUVQrd7dB6NGrbuoG586fxT9AIQQjxlBDiiBCi2s/nC4UQzUKIj42/+603U+OBq39hlj0O4SzdvsgqVY8n9ttni8ZyzAwZngEWBZF5T0o51fh7IHyzhi5Pvf8pTxlt1v0hgY8PnFAXYWsjdLbZbtNZYXnWSPUYZYdg5vxp/BPUIUgp3wUiX/2iAWD1nqOs3nM0oIxAcuB4+5m7coQTi542qQjhk0MtnkJZxqyH3RGLF2bOn8Y/VuUQ5gkhNgP1wN1Sym2+hIQQtwC3AJSWllqkenDx5M2zTMldMbnILUw/AHmVttp0xwVjPN9ITIfkHDge3QjB7PnT+MaKWYaNwEgp5RTg18Bf/QlKKZ+QUs6UUs7My4te4mlQ4MohgL1heqBZhOyRUY8QNOERtkOQUp6UUp4yni8H4oUQuWFbNkR54t09PPHuHv8CxgW5Yf8x1bosJj7iF6Evm17f5mP7tmgkOb0Iev40AQl7yCCEGA4cllJKIcRslJNpCvI1jR827jfXU6Ch+bRqwZ5VEvGL0NMm5ZB2N7byOW/BrFKoWWFbe3hfmD1/Gt8EdQhCiKXAQiBXCFEH/BiIB5BSPg5cA9wmhOgG2oFrpXRA65wByuM3zggsYJzay6cYKwwzI+8QPGzq66k4+mzBrJFqBWbrEbWjUxQIev40AQnqEKSU1wX5/BHgEcss0oRGVinsfD0Kin3s49g39VgbNYegCQ9duuwwfvPObn7zzu4AEuoO/eE+YyY4a6S6I3e122ST0v/3rT7WLbiSnFGcaQh+/jSB0A7BYWyvP8n2+pP+BYyQ/UjLafXafeoxQuyuO0xh9RPQ0933Xu0xHw6orxYheg4h6PnTBESvZXAYj1w/3ZTc5ZONHEK2W4Vg3lj/XwiDh2e1wNLfQu3lEJ8KwG0LfeQQElIhNS+qMw1mz5/GNzpCGHC48rXGGD6nQj0e2xtBlT2Gjk/P1u9NVmnUy5c1/Uc7BIfxqzd38as3dwWVW7PXKM9NzYOEtIg6hL6WZMfOzO+/stlP74Ps8sg6pyCYPX8a3+ghg8PY23gqsICRQzjW2qVeCwE5kb0IDzcb+YJje/v0H2ru8C2cUwHbXobu0xCXGDGb/BH0/GkCoh2Cw/jltdNMyV02ufDMi5wKOOxz+YglfH1BBbyIMWRQfOt8HzkEgGGjQPaqmYYI5TQCYfb8aXyjhwwDDh9j+JwKdQG6zQJEROexvepiD0TOqDOymgGHdggO4+E3anj4jRr/Aq59EXa7LfHNqYDeLjgZmSanr35stFfvaoNTag3Dnzf60dWX5IzOeoKg508TEO0QHEZ9cwf1/sbnbpzs6DnzIsIzDcfaOt1eqAu9qa3Lt3BKDiRlQlN0HEJu7QrG7f5dVHQPBnQOwWH875emBJFQEcKlk7xyCKAcwqgLLLfpK/NGgqvuybjQbzm3wrewEGrYEKUhw01Zm2HXSpA/P7tVvCYoOkIYDKQNh7hkj6RfxHDpCHSx5VREbciAlHD6JJzSnZf7g3YIDuN/VnzC/6z4xL+AkUN4Z6dbDiEmJqJTj3/d5JYvMC70Fz8KUCo9bJRqyd59OiL2BGLHIaNsuUnXIvQH7RAcxom2Tk64j9n90N7V4/lGTkXEHELraWP2ImMEtDQA0HK6x/8XctymHm2mq9uYBTmqHUJ/0DkEh/HfVwfbp1FFCJdUFXq+nVOuxs4RaE5yw5xSqEU5HWND128sKPf/Bfechs21CJOLM6EZaNIrHvuDjhAGGq7eM95j+GGjVXOSSHZgHuZWjBQohzDMVYsQhTyC6/xoh9AvtENwGP/19+3819+3B5V7c4dXT8Nc4058dKflNr288YCnDuC5dQFWNKbkQFJWVC7KHQ3GFnN6yNAvtENwGB1dvXR0BaoGVHfAs0RyjTbsjdYX5fSNy93C/87uIF3ycsdG5aLscZ2X4/ugO3guRuOJziE4jP/8fJUpuUUTvVqUpQ6DlGFw1HqH8OVZJbAflVQ0+FqgHAKofSJ2rrDclmBUFaWrHILsUU4hCuspBjI6Qhho+MshgIoSGq0fMpzRGcLPJa9SbTPXFsVNv/TUY8hoh+Aw/uNv23jg1S1B5d7YfujsN/PGqgjB4qbXy/rWLQhIV7Mbz6wJMqWYN049RmAIE4jtDSc5GW9sAqQTiyGjHYLDmNS0gru2fB66/K1nkMZ//UQI7cehzdptMQRuDsaYaYjvDVJ01JfktH+hUUdsmmocoxOLIaMdgsO4ukKS3t3kP9w17v6fm1h49meu8bLFd+UvTitWT4SAgokA3DAhIfCXMksgPsX2CGHC8DTyM5KN4ZNe9Rgq2iE4DuNufGRH6F91zTRE8q58wY9g7u1Q9cXAcjExkDsmehdl/nh1DvWeQSGhHYLDWLXdWJRzxF8tgvqBL9/qI4eQMUJ1RbY4sfjSelfNgYDEdP694wb+fbmJhVRRuEtvbzjJoZOnlUPobFFrKjSm0Q7BYcTFGk+CRAhxsT5yCDExkDva8gghPs74mRgzG0nxMSTFm/jp5FWqpi2nWyy1JxCxAoQQkD9BvdGfSGsIox2Cw1g41siQ+/shGyHwxb5yCKCy+xbflT8/pcjj9b9dNoF/u2xC8C/muYYwEZgK9UNlQToFGUmQb8xy+I20NL7QDsGpnNgPpwN1EPazliB/glqA1H48ImaFhGvq8UiA5dyWIwEBydmQXqQjhBDRDsFhrNzecOZFgDv937b42RdhuFHpaGEX5j+t96w5uPflLdz7cvBaCbLLITbR1rv09oaT1LvaxueP1xFCiGiH4DCS4mPPvPD1YzaGDEkJfqrOC6x3CCkuXUYOISslgayUINOOALFx6qI8tNUyW4IRHyuIdS3/zh+vhiu9AXo3aDzQaxkcxrmjc6EOiEuCRv+h9kXjC3x/kFag1jQcrrbMpssnFYLbSuYfLhpn/svDq6DmH8qR2dDjcExeKpw0NojJnwDdHWpNg2tJtiYgOkJwKnmVfsLdIHsruoqHDlnnEILqDETBJFU52eJjmjRiGHbmj1ePethgGu0QHMYb24wcQv7EgGH/XzYFmF8vmKSSaRaFyi98ZNQhGHf4u1/azN0vbTb35b6chpUOyj/b65upO2GUfedVAiKiu1oNNrRDcBjpicYornCy2hTF+85q5BAykgOM4QsmQne7ZV2Ys5LjPV4XZSZRlJlk7stGqbNdeYSk+BjiY42fdUKqGirYmMMY6OgcgsOYN2qYyiEUGvszNGyB9OFnyX3WXw4BzlyEh7eqQqUwWTRxOOwGVyj+LxdXmv9ycrZa12BThFCRmwqnWs+8UTgFDnxoi+7BgI4QnIar9t41W3DIOzQ3UZufNw5ErIWhcpjrAQqqLM5phEDhFNVnstXaFaCDlaAOQQjxlBDiiBDC57+oUPxKCLFbCLFFCDHdejOHDq9vO0QvApIyVDvzBi+HYDiMlzYc9H+Q+CS1sMiiUHnpOqMOwcghfO+FTXzvhU3mDzC8Sq3e7Gq3xJ5AbK9vpva4mx5XpHWWY9X4wkyE8AywKMDnlwBjjL9bgMfCN2vokp0SfyaXXzjlbIdgkJuWGPhAhVOh/mNLbMpLd+lSllXkpVGRl2b+AAVVap8GG6oGUxJiSYhzq+UYbrS1bzBRSKUJ7hCklO8CgfpgXQU8KxVrgSwhhJ9Ce00wZpfnqMU5oBKLJ2q9ypBVhPCZcfmBD1Q0DU4dgpMNgeVMcOF4T13f+ewYvvPZMeYPUDRVPdaHEFX0k7JhKQzPcEt4puRAVqlfx6rxxIocwgjObAUKKiU2wo+sJhju6/fdE4veBCvyKZqmHus3WmOXGZ3+yBoJyTnW2uIPXwVQASItjSdWOARfvxKfWSghxC1CiPVCiPWNjY0WqB58rNh2iB7X2RvucghuP2bDYSz9MMC+CADDJ6nEogV35T+u81zLcMfzG7nj+RAubiFgxHQ4GPkIYXvDSfY1tXq+WThFbRrTcTLi+gc6VjiEOqDE7XUx4HPljZTyCSnlTCnlzLy8PAtUDz7y0hLO3OFSh0FGsdfdTTmEwqyUwAdKSFGVegfDvyufqTlQdk0oymBCUUaIB5kOjTugszW4bBikJcaSFO81m+5yrId0HiEYVjiEV4GbjNmGuUCzlDL8gesQZUZpFrHuIe+IaXBw/VlyfX0TAlE0TUUIYbYR+0ylp67bF47m9oUh1jeMmKESixEO3Uuzkz1zCHBm+GSBcxzsmJl2XAqsASqFEHVCiK8LIW4VQtxqiCwH9qJKV34H3B4xa4cMbg6heLZanHPKGGIF2pfBm6Jp0H5MJSbDIRSd/hhhzEbbcVF625mWp/IYdR9FXvcAJ2ilopTyuiCfS+Dbllk0xFlR3cCFvfLMP0zxLPV4cD1UXtIn99za/SyeEeRgrouwfiNkj+y3Tc+tq2Ux4HJUt/5hAwCP3xjMADfS8tXwJ8KJxR0NJ0nobeestY0ls2Hf+7atuhyo6EpFh1GQmXRm2hFUQiwmzu3upu7WJcNSgx8sfyLEJsDBDWHZNDLHMwSfPjKL6SOzQj/QiOlh2xKM9KQ4UhN93OeKZ0FLQ9929hrf6LUMDmNacSbUuTmEhBRV2ONVj3/+2CB1CABxCWrYULsuLJvOHZMHNfTdWW85r5+9BUZMhx2vqjLi1GFh2eSP4qwk6Pb1gRFpHfgQMosjonswoCMER+IV0pbMVmPv3p7QE4Qlc6Dh4wA7QZnAqr0NSuaqxwNrrTleKBRUqaYzdWcnaDVn0A7BYfyjuoHOXq8LsHgWdLUapb/qs6dX7zN3wNJ50NMZVj3Cc2tdupSj+saSj/jGkn4k6IqmqSFM7Zp+2xKMHQ0n2XXER3NaV7RUp1c+BkI7BIcxIiuJGO+kV/FM9ej2Yx6dl27ugCVz1GMYF6H3uoX5o3KZPyo39APFJ6l6hNrIRQiZyXGk++sVUTxTTXt2B9mXcgijHYLDmDwik7gYr3+W7HJIzVcXkhG+n2umDgHUWH3YGDjQ/zzC/Ioc9cRwVF9bUM7XFpT372Clc1W00tnWb3sCUZSZdHYdgouSOWFHS4Md7RAch4/xuhBQdo6aNutPf8PSucqZ9PaGaZsF03Uj50Nvd4RnG/zYWTpfPe57P4K6BzbaITiM5VsPcbrHh1MoW6CmzI7vA+B37+01f9DSudBxot87KD27Zp/H65uf+pCbn+rnWLxktnqM0LBhR8NJag772ToudZiaitUOwS962tFhjByWTKyvaHrkAvVo/JgnFGWaP2jpPPW4/4MzW5yFQOXwdDhO35DBezl0SCRnq/botav7f4wAqH4SAWZFyhbApj9Ad6dKNGo80BGCw5hYmEFcTOzZH+RVQkoufPoeAOeMDiGpl1Ohdob+9J/9smlOWY7H6xvnlXHjvLJ+HQtQw4badeqitJjhGUlqb0d/lJ0DXW06j+AH7RCciK/SWiHU3a1vDj+E8bwQUH6+cib9yiOEsS+DL8rPV9OoEckjGHs7+mPkOepx33sR0D3w0Q7BYSzfWk97l5/9FMoW9D197J97fMv4o+J8tdDpcOh9Fpes2aeeGI7qhifXcsOTYeQAys8FEQN73+7/MfzwSUMz2w8F6HuQmquGLDqP4BPtEBxGRW4qcbF+/lnKzu17OrU0O7QDl5+vHve+E7JNVV69Dy6fXMTlk4v8SJsgOVsVCfXDlmDkpCYwLDVIv8myBWoaVtcjnIV2CA5j3PB04r3rEFzkndkPYV5FiIVBGYWQWwl7Q88jzPByPtfNLuW62aUhH8eDioWqjNjiLkb56YmBcwgAoy5QeYQIVkwOVLRDcCL+lucKcWY9gOxHLqDifHURhJzMs2gtg4ctnwHZE53QvexcVUK9e5X9uh2OdggOY/mWeto6fXmoT4gAABawSURBVC3XMxh1AQCvrXwj9INXLOzXnfEZ17oJw1F9+bdr+PJvw7y7lsyGuGTL8wifNDSztd5PHYKLxDQ1FbtLOwRvtENwGKML0oiL9THt6GLW1zmWWUX8jBtDP3j5+erOuCs0ZzK12LPm4ZoZxVwzI8wlxHGJaiy/e5V1qylR+1XkBxsyAIy+UPV4bA6wae4QRDsEhzE2P42EuAD/LKm55Nz1AZ+74ILQD56YBuXnQc0/Qvra1BKXQ1ARwpdmlvClmSX+v2CWykVwbC8c3RX+sQxy0xIoSA+SVAQYc5F61MMGD7RDcBg9UgYdsXf19NLV0891CWMXqZbkR3ebt6nXs6diWPq9bQGoWR7+sQx6Za/aCi8YeeNUsdaulZbpHgxoh+Aw3qhuoPW0nzoEg8VPrmPxk/1cvTjmYvW4c4Xpr/yxrx+CBfrdySxWW62FYEswdh0+RfXB5uCCQqhzseft8JrHDDK0Q3AYYwuCDBmAa2eXcO3sfobs2SNVYU4IF+G0Ulf/RBG+fm8qL1E1ARbtzpyXlkBBZrI54fFXqIrJCBRIDVS0Q3AYo/LSPDcr9cEXphXzhWlhJPXGLoL9q6Et0JadZ5g0wjOpGLZ+dyovUVOoISY6/ZGTajKHAGr6MSkTtr9qie7BgHYIDqO7pzdoDqG9s4f2zsDDioBMuErVAHzyminxvnyBkUMIW787hVPVWH6HNRdlb28vvlaP+yQuAcZeonIYPV2W6B/oaIfgMFbtOExLR+CL7StPf8hXng6jN2DhFNWFqfplU+IveO0jGbZ+d4SAiV9Q2f72E2Efbk9jK1vrQ6h+HH+F6hWh1zYA2iE4jsqCNBLjAw8ZFs8dyeK5/d94pe8i/PRdaD0aVHyGVw4hbP3eTLxatTb75O9hHyo3LYEiszkEUIVe8Smw/ZWwdQ8GtENwGOW5qSQGSSpeMaWIK6aEsbgIoOpqNWwwEar3bexqDBks0e/OiOlqq7XqP4d9qOzkOHOFSS4SUlQeY/tfI9KfYaChHYLD6OzuCTqPfrKji5MdYY55C6pg2GhTw4aOLs9Sakv0uyMEVH1RrX40EbEEokdKur3b2Adj8rXQftyyxOZARjsEh/FOzRFaglxs31yynm8uCXPDESFg0pfU2Pn4/oCif97g2v5MWKffm6ovqohl21/COsy+o6fYejDEFZSjLoDUPNjyQli6BwPaITiMccPTSYoP3Oryq+eU8dVzysJXNvV69bh5aUCxmV77OFqm353hVTB8Emx8NqzD5KYlMiI7JbQvxcZB1TWw83UVKQxhtENwGKXZyUFzCIuqCllUVRi+sqxStSR60x8DtlarLDA2ajFyCJbp92b6zXBoC9R/3O9DZCbFkW+2DsGdKV9WiU2TMy+DFe0QHEZHdy/BhsDHWjs51mpRAmzajdBcC/ve9SvS1tfSTViv351J10BsouqK3E+6e3vpNF2I4EbhVJVX2fC0pasvBxraITiM93cdobkjQD8E4LbnNnDbcxY1KB13marW27DEr8grmzyXCFuq353kbJhwJWx5Cbra+3WI2qY2tppZy+CNEDDr63BoK9T1Y9/KQYJ2CA5j3PAMUhIC5xC+eW4F3zy3whqF8ckw9QY1/Xiy3qfIrDLPrdws1e/N9JvhdDNsXdavr+emJVCak9o/3ZP+HySkw0e/79/3BwHaITiM4qwkEoOsZbhwQgEXTiiwTunsW9RW834uhNF5nheY5frdKVugQve1j/UrdM9IiiOvPzkEUP0iplwL2162bLHVQEM7BIfR3tVDT5AL4UhLB0daLFyym1OuinM2PO1zKfCpvpZuIjL63REC5t4OR7b1a2OZrp5eTneH0ath1tdVcnH90IwStENwGKt3H6W5PXAO4c7nN3Hn8xbvPDTnVmhrgq1/Ouuj1z72HEpERL87VV9UdQFrfhPyV+uOt7Ml1DoEd/LHw5jPqQils7X/xxmgaIfgMCYUppOaGDiHcNvCUdy2cJS1isvPU81K3v8F9Hg6pDnlRht2I4cQEf3uxCfBrG/Crtfh8LaQvpqXFk95bj9zCC4W3KU2tdn0XHjHGYCYcghCiEVCiBohxG4hxD0+Pv+KEKJRCPGx8fcN600dGhRmBs8hLKzMZ2FlGBuu+kIIOP+HqsdhtWdCrzw3zSUUOf3ezP4mJGbAP/8npK+lJcaRmxbCWgZfjJynujKv/vWQWxYd1CEIIWKBR4FLgAnAdUKICT5EX5RSTjX+nrTYziFDa2d30PX89SfaqT/Rv2m5gFReqhJ67/5cJRkNTnZ41hxETL87KTlqGLP9FThUbfprXT29/rfCC4UFd0HzgaBVnIMNMxHCbGC3lHKvlLITeAG4KrJmDV3W7W3ieFvgop+7XvyYu17sfzWfX2Ji4PwfQNNuj5WH/9jaoJ4YQ4aI6fdm3u0hRwn1J9pCX8vgizEXw4iZ8PZ/97smYiBixiGMAA64va4z3vPmi0KILUKIZUIInw33hBC3CCHWCyHWNzY29sPcwc/EwgzSkuIDytx5wRjuvGBMZAwYdwUUTII3/7NvxmFuxTD79LuTnA1zb1M1EnXmFlPlpiVSkZ8WXDAYQsBF/wEt9bDut+Efb4BgxiH4WovrHdT+DSiTUk4GVgE+y96klE9IKWdKKWfm5eWFZukQoSAjiaQgOYQFY3JZMCbEvR3NEhMDix5S5cxrHwVgZI5rsZCIvH5v5t8JaQXwjx+a2so+NSGW3NQEa3SXLVAzDu89bLr/5EDHjEOoA9zv+MWAxzyUlLJJSunaSvd3wAxrzBt6nDrdRXeQOoTapjZqm9oiZ0T5eVB5mboQWg5zwmsIE3H97iSmw2d/DAfXw9aXgop39vTS1mXBnhEuLvwxdLbAWw9ad0wHY8YhfASMEUKUCyESgGsBjzY7Qgj3pW9XAjusM3FosX5fE8fbAme2v79sM99ftjmyhlz8n2q79JX/zsrtnjkEW/S7M+U6tX38qh8H3S36cHN7/9Yy+KNgIsz+Fqx/yvSwZSAT1CFIKbuBO4DXURf6n6SU24QQDwghrjTEviOE2CaE2Ax8B/hKpAwe7EwsyiQjSA7hrovGctdFYyNryLBRKtO+5UUuSdyq3nMlFe3Q705MDFz6v9ByCFbeH1B0WGoCo/PTrdV/wb9BeiG89r2zajQGG4ErYAyklMuB5V7v3e/2/F7gXmtNG5rkpSVAkCar3km+iHHe3bDjVdIat0RHvzvFM2Het2HNIzDx82onax+kJMSSYlUOwUViOlzyU/jTTbD6/+Dcf7X2+A5CVyo6jJaOLrqCFCLsaTzFnsZTkTcmLhGufCR6+r254EeqD+Srd/odOnR293DKqj0j3Bl/pepU/fZDcDACS78dgnYIDmNT7XGOBckh3PfyVu57eas9BpXM4s3kRRyNOTOrYKt+d+KT4fOPQfNBePUOn6shj7SctqYOwRsh4PJfqBmPP38TTkfBIdqAqSGDxj4mFmWQdjhwDuEHiyptskaRde1j7Je9uFyC3fo9KJmtMv8r71f1AXNv9fh4WGo8KdkW5xBcJGfD1U/AM5erKOWap/ryKoMF7RAcxrDU4DmEGSNzbLLGtz679Z/FvDth/xp440eqOWvZgr6PkuNjSU7tZz8EM5QtUA5p1U/UDMR5d0dOVxTQQwaHcbK9M2hPwJpDLdQcarHJorP12a3/LGJi4AuPQU4FvHA9NNb0fXS6u4eTQVrQhc0531Mt7N96cNDt+KQdgsPYfOBE0Aam979Szf2vmF/wEy7e+uzW75PkbLjhJdWU9blroFn1fWxq6aA6EjkEd4SAK38NxbPgz9+A3W9GVp+N6CGDw6gakUFKY+CQ975Lx9tkjW99duv3S/ZIuOFPsORKeOYyuPk1ctISSR6WGfy74RKfrHQ/cwW8cAPc+DKMnB95vRFGRwgOIzslIei+DFNKsphSkhVQxkq89dmtPyBF0+DGv0LbcXjmUpJaaslOCZyUtYzkbLjxL5BZDH+4Wm30MsDRDsFhnGjr5HSQHMK2+ma21VtYnhsEb3126w9K8Qy46S/Q0Qy93bQ21QX/jlWk5cFXl0PeWFh6ndr0ZgCjHYLDqD7YTNOpwDmEB/62nQf+tt0mi87WZ7d+U4yYAV9fSRfxrDpZaq/utHz4yt+h/Fx45Xa1MnOA7iStcwgOY9KIDJKOBs4h3H+Fr4ZVkcNbn936TZNXyc5v7GS0iMJ9LjEdblgGK3+slo3Xb1I1C9ll9tsSBtohOIzM5HgIkkOYWGRD0iyAPrv1h8LE4ijWSMTGq14SxTPgb9+D38xXTVZmfl1NlQ4ABoaVQ4jjrZ10dAfOIWw+cILNB07YZNHZ+uzWHwqOsK3qi3D7GiidC8vvht9fBLXromuTSbRDcBg7GoLnEB5avoOHltvXcsJbn936Q8ExtmUWw+I/w+cfVzUST12sVkuG0DA2Gughg8OoGpFJYlPgHMIDV1XZZI1vfXbrDwVH2SYETL1ObWC7+hH44P9UZePoi+Cc70DZuY5bC6EdgsPISIyFID0VK4dHaPGOSX126w8FR9qWkAoLf6j2mlj/e1j7OCy5ArLL1Ua7U66FLJ99iW1HDxkcxrG2zqD7CmzYf4wN++1r+umtz279oeBk20jJgfO+D3dVq6FEZjG8/SD8sgqe+Az88+dqSNGPTW6tQkcIDmPnoZOM6O4k0P3iZyvUYp4XvzXPFpu89dmtPxScbFsf8clqKDH1Ojj2qdoDo2a5cg5vPwgpuaoMeuR8lZjMn6Ca1diAkFHyRjNnzpTr1w/+ppWhcuoP15NwbBcJ3/3Ir4yrW9GoPAv2HzCBtz679YeCk20LSssh2PUG7F8N+z+AE7Xq/Zg4yK1US70LJqpVntnlqsYhMfT/TyHEBinlTF+f6QjBYaQlBM8h2P1j99bn5IvNybYFJX04TL9J/QGcOAB1H8HhajWU+PQ92PKi53dS89Uir4wiyBihWtbH939vS2c7hNOnVH36EOL4yWYSu3pICSCzdm8TYF+zU299dusPBSfbFjJZJeqv6uoz77WfgOOfqk15j32qnp+ohcPblcP43ENhqXS2Q9jyIvz9X6Jtha1kA3vjRlMRQOYXK3cC9o2TvfXZrT8UnGybJSRnQfI0tcozAjg7h9BYA7Vr7THIITS1dtKZP4XCcbP9yrh2TSodFiiOsA5vfXbrDwUn2+YUBm4OIa9S/Q0hzAS6dv/YvfU5+WJzsm0DAV2H4DDe33WU93cdDVvGSrz12a0/FJxs20DA2RHCEOTXb+0CCLi7shmZSNpkt/5QcLJtAwFn5xCGIPUn2gEoykoOSyaSNtmtPxScbJtTGLg5hCGImR+y3T92b31OvticbNtAQOcQHMY7NUd4p+ZI2DJW4q3Pbv2h4GTbBgI6QnAYj72zB4CFlflhyUTSJrv1h4KTbRsI6ByCwzjS0gFAfrr/8lMzMpG0yW79oeBk25yCziEMIMz8kO3+sXvrc/LF5mTbBgI6h+AwVm0/zKrth8OWsRJvfXbrDwUn2zYQ0BGCw/jde3sBuHBCQVgykbTJbv2h4GTbBgI6h+AwXBu95qQmhCUTSZvs1h8KTrbNKegcwgDCzA/Z7h+7tz4nX2xOtm0goHMIDmNFdQMrqhvClrESb3126w8FJ9s2EDAVIQghFgH/B8QCT0opf+r1eSLwLDADaAK+LKXcZ62pg5vH/7mHycWZPP3BPgAWVRWyes9RttQ1c+v5o0zLRNKmfU1txMbAs2v2U5SZHHH94djqJNsGFFLKgH8oJ7AHqAASgM3ABC+Z24HHjefXAi8GO+6MGTOk5gwf7G6U0x54Q67c3iCb2zv7Xn+wuzEkmUja9MS7u2XZD1+Tv3qzxhb94djqJNucBrBe+rkugyYVhRDzgJ9IKT9nvL7XcCT/7SbzuiGzRggRBxwC8mSAg+uk4tms3nOUO57fxOI5pTy3rpZHrp/G/FG5IctE0qbbFlbw2Dt7bdMfCnafm4FKoKSimRzCCOCA2+s64z2fMlLKbqAZH70+hBC3CCHWCyHWNzY2mrF9SDF/VC6zy7L51Vu7WTyn1OeP2YxMJG365rmjbNUfCnafm8GIGYfga68p7zu/GRmklE9IKWdKKWfm5eWZsW9IsXrPUVbtOMKIrCSeW1fL6j1nN/owIxNJm3733h5b9YeC3edmUOJvLCHP5AfmAa+7vb4XuNdL5nVgnvE8DjiKUePg70/nEDxxjXnf+uSQbDvdHTCHEEgmkja5cgiPvr3TFv3h2Ook25wGAXIIZiKEj4AxQohyIUQCKmn4qpfMq8DNxvNrgLcMxRqTbKlr5pHrp/GZygKSE2KZPyqXR66fxpa65pBkImlTTy/cd9k4YkSMLfrDsdVJtg0kTFUqCiEuBX6JmnF4Skr5X0KIB1Ce5lUhRBLwB2AacAy4Vkq5N9AxdVLRN3/ZVAfAF6YVhyUTSZvs1h8KTrbNKYRdqSilXA4s93rvfrfnHcCXwjFSo3jhQ5W/DfSDNiMTSZvs1h8KTrZtIKDXMjiMrp5eAOJj/Y/mzMhE0ia79YeCk21zCnotwwDCzA/Z7h+7tz4nX2xOtm0goM+ew3hp/QFeWn8gbBkr8dZnt/5QcLJtAwHtEBzGsg11LNtQF7aMlXjrs1t/KDjZtoFA1HIIQohGYL8J0VxUXYOT0TaGj9PtA+fbaNa+kVJKn5WBUXMIZhFCrPeXAHEK2sbwcbp94HwbrbBPDxk0Gk0f2iFoNJo+BoJDeCLaBphA2xg+TrcPnG9j2PY5Poeg0WjsYyBECBqNxia0Q9BoNH042iEIIRYJIWqEELuFEPdE2x53hBAlQoi3hRA7hBDbhBDfjbZN/hBCxAohNgkhXou2Lb4QQmQJIZYJIT4xzue8aNvkjhDiLuPfuFoIsdRY3Rttm54SQhwRQlS7vZcjhFgphNhlPGaHelzHOgQhRCzwKHAJMAG4TggxIbpWedAN/KuUcjwwF/i2w+xz57vAjmgbEYD/A1ZIKccBU3CQrUKIEcB3gJlSyipUC4Bro2sVAM8Ai7zeuwd4U0o5BnjTeB0SjnUIwGxgt5Ryr5SyE3gBuCrKNvUhpWyQUm40nregfsTevSajjhCiGLgMeDLatvhCCJEBnAf8HkBK2SmlPBFdq84iDkg2GginAPVRtgcp5buo3iPuXAUsMZ4vAT4f6nGd7BDMNHd1BEKIMlRzmHXRtcQnvwR+APRG2xA/VACNwNPGsOZJIURqtI1yIaU8CPwvUAs0AM1Syjeia5VfCqSUDaBuWEB+qAdwskMw1bg12ggh0oA/A9+TUp6Mtj3uCCEuB45IKTdE25YAxAHTgceklNOAVvoR6kYKYxx+FVAOFAGpQojF0bUqcjjZIdQBJW6vi3FAqOaOECIe5Qz+KKV8Odr2+OAc4EohxD7UkOsCIcRz0TXpLOqAOimlK7pahnIQTuFC4FMpZaOUsgt4GZgfZZv8cVgIUQhgPB4J9QBOdghmmrtGDSGEQI17d0gpH462Pb6QUt4rpSyWUpahzt9bUkpH3d2klIeAA0KISuOtzwLbo2iSN7XAXCFEivFv/lkclPT0wr3Z8c3AK6EewLEdk6SU3UKIO1At3l3NXbdF2Sx3zgFuBLYKIT423rvP6D+pCY07gT8ajn8v8NUo29OHlHKdEGIZsBE1s7QJB5QwCyGWAguBXCFEHfBj4KfAn4QQX0c5spD7nOrSZY1G04eThwwajcZmtEPQaDR9aIeg0Wj60A5Bo9H0oR2CRqPpQzsEjUbTh3YIGo2mj/8PLHFVIVX0PMsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "background = 0.2\n",
    "kernel=lambda dx: 2**(-dx*2)\n",
    "\n",
    "x = Hawkes(np.random.RandomState(0), background, kernel, 0, 10)\n",
    "lamb = lambda xnew: background + np.sum([\n",
    "    kernel(xnew-xi) if xnew>=xi else 0 for xi in x])\n",
    "\n",
    "pl.plot(x, np.zeros_like(x), 'x', label='sample')\n",
    "[pl.axvline(xi, ls=':') for xi in x]\n",
    "pl.plot(np.linspace(0,10,1000),\n",
    "        [lamb(xnew) for xnew in np.linspace(0,10,1000)], label='$\\lambda$')\n",
    "pl.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Environment (conda_mxnet_p36)",
   "language": "python",
   "name": "conda_mxnet_p36"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
